# 所有的数据处理函数
# 张浩 2019年3月24日

import struct
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from datetime import datetime
from scipy.interpolate import interp1d
from scipy import interpolate
import seaborn as sns

def filter_Remove_signal_distortion(signalIn, thresh, searchLen):
    """  去除由于相位展开而造成的信号畸变

    搜索signal中的相邻两个采样值的差，如果有一个正差大于门限thresh 则表示为发现波形畸变
    标志位置1  后向搜索至searchLen长度  如果发现有负跳变  则表示找到波形畸变结束  结束标志置1
    将开始和结束之间的数值整体减掉跳变值  返回修改后的信号

    :param signalIn:待处理信号
    :param thresh:门限
    :param searchLen:向后搜索的长度
    :return:处理后的信号signal np一维数组
    """
    startIndex = 0
    endIndex = 0
    i = 0
    Flag = 0
    oldValue = signalIn[0]

    signalInLen = len(signalIn)
    while i < signalInLen:
        diff = signalIn[i] - oldValue
        if diff > thresh:
            startIndex = i
            xiuzhengVakue = diff
            searchLenRel = min(searchLen,(signalInLen-i-1))
            ii = 0
            while ii < searchLenRel:
                diff = signalIn[ii + startIndex] - oldValue
                if diff < -thresh:
                    endIndex = ii + startIndex
                    Flag = 1
                    break
                oldValue = signalIn[ii + startIndex]
                ii = ii + 1
            if Flag:
                jj = 0
                while jj < endIndex - startIndex:
                    signalIn[jj + startIndex] = signalIn[jj + startIndex] - xiuzhengVakue
                    jj = jj +1
                Flag = 0
                i = endIndex-1
            else:
                jj = 0
                while jj < signalInLen-startIndex:
                    signalIn[jj + startIndex] = signalIn[jj + startIndex] - xiuzhengVakue
                    jj = jj + 1
        oldValue = signalIn[i]
        i = i+1

    i = 0
    while i < signalInLen:
        diff = signalIn[i] - oldValue
        if diff < -thresh:
            startIndex = i
            # print('负跳变：%d' % startIndex)
            xiuzhengVakue = abs(diff)
            searchLenRel = min(searchLen, (signalInLen - i - 1))
            ii = 0
            while ii < searchLenRel:
                diff = signalIn[ii + startIndex] - oldValue
                if diff > thresh:
                    endIndex = ii + startIndex
                    Flag = 1
                    break
                oldValue = signalIn[ii + startIndex]
                ii = ii + 1
            if Flag:
                jj = 0
                while jj < endIndex - startIndex:
                    signalIn[jj + startIndex] = signalIn[jj + startIndex] + xiuzhengVakue
                    jj = jj + 1
                Flag = 0
                i = endIndex - 1
            else:
                jj = 0
                while jj < startIndex:
                    signalIn[jj] = signalIn[jj] - xiuzhengVakue
                    jj = jj + 1
        oldValue = signalIn[i]
        i = i + 1

    return signalIn

def getDataFromFile(fileName,fs=200,len=1800,time=[0,120]):
    """  从文件中获取胸部位移数据

    :param fileName: 数据文件的文件名
    :param fs: 信号的采样频率 Hz
    :param len: 信号的总长度 秒
    :param time: 要获取的信号的区间
    :return: np数组  胸部位移的值
    """
    file = open(fileName,mode='r')
    fulldata = np.zeros(len*fs)
    distance = np.zeros((time[1]-time[0])*fs)
    i = 0
    while True:
        temp = file.readline()
        if temp:
            fulldata[i] = float(temp)
            i=i+1
        else:
            break

    i = 0
    while i < (time[1] - time[0])*fs:
        distance[i] = fulldata[time[0]*fs+i]
        i = i + 1
    file.close()

    return distance

def parseFileToTXT(originalFile,parsedFile,fs=200,len=1800,time=[0,len],packageLen = 160):
    """  将板子上传来的bin文件转换为txt类型的列向量

    :param originalFile: 原始文件 .bin
    :param parsedFile:   转换后的文件
    :param fs:   采样频率
    :param len:   信号总时长
    :param time:  要转换的区间长度  默认全部转换
    :return:
    """
    time1 = datetime.now()
    dataOut = np.zeros((time[1]-time[0])*fs)
    dataAll = np.zeros(len*fs)
    package = []
    with open(originalFile,mode='rb') as ofile:
        while True:
            temp = ofile.read(packageLen)
            if temp:
                package.append(temp)
            else:
                break
    print('读取bin文件耗时：')
    time2 = datetime.now()
    print(time2-time1)
    i = 0
    for ii in package:
        temp = ii[64:68]
        temp_float = struct.unpack('f',temp)
        dataAll[i] = temp_float[0]
        i = i + 1
    print('将bin文件转换为float耗时：')
    time3 = datetime.now()
    print(time3-time2)

    dataOut = dataAll[time[0]*fs:time[1]*fs]

    print('切片耗时：')
    time4 = datetime.now()
    print(time4 - time3)

    if parsedFile:
        with open(parsedFile,mode='w') as pfile:
            i = 0
            while i < (time[1]-time[0])*fs:
                pfile.write(str(dataOut[i]))
                if i == (time[1]-time[0])*fs-1:
                    break
                else:
                    pfile.write('\n')
                i = i+1
    print('将数组写入TXT耗时：')
    time5 = datetime.now()
    print(time5 - time4)
    print('总耗时：')
    print(time5-time1)
    return dataOut

def filter_discha(signalIn):
    """
    计算连续的相位差
    :param signalIn: 输入的待处理的信号
    :return: 输入信号的连续相位差
    """
    discha = np.zeros(len(signalIn))
    i = 1
    while i < len(signalIn):
        discha[i] = signalIn[i] - signalIn[i-1]
        i = i + 1
    return discha

def filter_removeInpluseNoise(signalIn,thresh = 0.1):
    """
    去除脉冲噪声
    :param signalIn: 输入的待处理信号
    :param thresh: 门限
    :return: 去除脉冲噪声后的信号
    """
    i = 1
    while i < len(signalIn) - 1:
        forwardDiff = signalIn[i] - signalIn[i - 1]
        backwardDiff = signalIn[i] - signalIn[i + 1]
        if ((forwardDiff > thresh) and (backwardDiff > thresh)) or (
            (forwardDiff < -thresh) and (backwardDiff < -thresh)):
            signalIn[i] = signalIn[i - 1] + (0.5 * (signalIn[i + 1] - signalIn[i - 1]))
        else:
            signalIn[i] = signalIn[i]
        i = i + 1
    return signalIn

def filter_weifen2(signalIn):
    """
    仿照文献Emotion Recognition using Wireless Signals中的方法
    将信号做二次微分 求解信号的加速度  抑制呼吸  增强心率信号
    :param signalIn: 输入的待处理信号
    :return: 二次微分后的信号
    """
    disDaoshu = np.zeros(len(signalIn))
    discha = signalIn
    i = 3
    while i < len(discha)-3:
        disDaoshu[i] = -(
        4 * discha[i] + (discha[i - 1] + discha[i + 1]) + 2 * (discha[i - 2] + discha[i + 2]) + (discha[i - 3]) +
        discha[i + 3]) / (16 * 0.004 * 0.004)
        i = i + 1
    return disDaoshu

def average2D(data):
    """
    计算平均模板  即二维数组的每列的平均值
    :param data: 输入的二维矩阵
    :return:
    """

    sigLen = len(data[0])
    dataOut = np.zeros(sigLen)
    i = 0
    sum1 = 0
    while i < sigLen:
        ii = 0
        while ii < len(data):
            sum1 = sum1 + data[ii][i]
            ii = ii + 1
        dataOut[i] = sum1/len(data)
        i = i + 1
    return dataOut


def computerTemple(signalIn,fs):
    """
    计算心跳信号的模板
    :param signalIn:输入的样本信号
    :param fs: 采样率
    :return: 信号模板 0.8秒 采样率200Hz时为160采样点长度
    """
    signalLen = len(signalIn)
    menxian  = 0.75
    #1、寻找出所有的极大值点
    peaks,_  = signal.find_peaks(signalIn)

    #2、截取极值点前后0.4秒的数据作为心拍  放入二维数组
    beatLen = int(fs*0.4) # 80点
    # 计算beats 的长度   看看有多少个峰值点满足要求
    i = 0
    for beatIndex in peaks:
        if beatIndex < beatLen or beatIndex > signalLen-beatLen:
            continue
        else:
            i = i + 1
    beats = np.zeros([i,beatLen*2])
    i = 0
    for beatIndex in peaks:
        if beatIndex < beatLen or beatIndex > signalLen-beatLen:
            continue
        else:
            ii = 0
            while ii < 2*beatLen:
                beats[i][ii] = signalIn[beatIndex-beatLen+ii]
                ii = ii + 1
            #print('beats[%d]:' % i)
            #print(beats[i])
            i = i + 1
    #print('len(beats):%d' % len(beats))
    # 3、计算每两段信号的相关性  得到相关性矩阵
    corr = np.corrcoef(beats)
    #print('相关系数矩阵：')
    #print(corr)

    # 4、 根据相关性矩阵做聚类分析
    m = 1
    julei = [[0],]
    #print('开始聚类算法,beatsNum:%d' % len(beats))
    for i in range(1,len(beats)):
        #print('No.%d' % i)
        #print('此刻聚类的类个数：%d' % len(julei))
        #print(julei)
        m = 0
        minInC = []
        while m < len(julei):
            # 在每个已存在的聚类中计算beats[i]与该聚类样本间的相关性最小值
            temps = []
            for mm in julei[m]:
                temps.append(corr[i][mm])
            minInC.append(min(temps))
            m = m + 1
        #print('minInC:')
        #print(minInC)
        # 选出最相关的一个聚类  得到相关性 和聚类的索引
        maxXiangguanJL = max(minInC)
        maxXiangguanIndex = 0
        ii = 0
        while ii < len(julei):
            if minInC[ii] == maxXiangguanJL:
                maxXiangguanIndex = ii
                break
            ii = ii + 1
        #print('最相关的一个聚类：%d' % maxXiangguanIndex)
        #print('最大相似度：%f' % maxXiangguanJL)
        if maxXiangguanJL > menxian:
            # 属于已经存在的聚类  加入到该聚类
            julei[maxXiangguanIndex].append(int(i))
        else:
            # 新生成一个聚类
            temp = [i,]
            julei.append(temp)
    #print('聚类结果：')
    #print(julei)

    juleiLens = np.zeros(len(julei),dtype=np.int)
    i = 0
    while i < len(julei):
        juleiLens[i] = len(julei[i])
        i = i + 1
    maxNum = max(juleiLens)
    i = 0
    while i < len(juleiLens):
        if maxNum == juleiLens[i]:
            maxIndex = i
            break
        else:
            i = i + 1
    dataTemp = np.zeros([juleiLens[maxIndex],beatLen*2])
    i = 0
    while i < juleiLens[maxIndex]:
        dataTemp[i] = beats[julei[maxIndex][i]]
        i = i + 1
    #sigModel = average2D(dataTemp)
    sigModel = dataTemp[0]
    #print('maxIndex：%d signalNum：%d' % (maxIndex,len(julei[maxIndex])))
    #print(sigModel)

    return julei,beats,sigModel

def computerTempleJpeck(signalIn,fs):
    """
    计算心跳信号的模板  以J峰为开始
    :param signalIn:输入的样本信号
    :param fs: 采样率
    :return: 信号模板 0.8秒 采样率200Hz时为160采样点长度
    """
    signalLen = len(signalIn)
    menxian  = 0.75
    #1、寻找出所有的极大值点
    peaks,_  = signal.find_peaks(signalIn)

    #2、截取极值点前后各0.4秒的数据作为心拍  放入二维数组
    beatLen = int(fs*0.4) # 200点
    # 计算beats 的长度   看看有多少个峰值点满足要求
    i = 0
    for beatIndex in peaks:
        if beatIndex < beatLen or beatIndex > signalLen-beatLen:
            continue
        else:
            i = i + 1
    beats = np.zeros([i,beatLen*2])
    i = 0
    for beatIndex in peaks:
        if beatIndex < beatLen or beatIndex > signalLen-beatLen:
            continue
        else:
            ii = 0
            while ii < 2*beatLen:
                beats[i][ii] = signalIn[beatIndex-beatLen+ii]
                ii = ii + 1
            #print('beats[%d]:' % i)
            #print(beats[i])
            i = i + 1
    #print('len(beats):%d' % len(beats))
    # 3、计算每两段信号的相关性  得到相关性矩阵
    corr = np.corrcoef(beats)
    #print('相关系数矩阵：')
    #print(corr)
    # 热力图
    # plt.subplots(figsize=(9, 9))
    # sns.heatmap(corr, annot=False, vmax=1, square=True, cmap="Blues")
    # plt.show()

    # 4、 根据相关性矩阵做聚类分析
    m = 1
    julei = [[0],]
    #print('开始聚类算法,beatsNum:%d' % len(beats))
    for i in range(1,len(beats)):
        #print('No.%d' % i)
        #print('此刻聚类的类个数：%d' % len(julei))
        #print(julei)
        m = 0
        minInC = []
        while m < len(julei):
            # 在每个已存在的聚类中计算beats[i]与该聚类样本间的相关性最小值
            temps = []
            for mm in julei[m]:
                temps.append(corr[i][mm])
            minInC.append(min(temps))
            m = m + 1
        #print('minInC:')
        #print(minInC)
        # 选出最相关的一个聚类  得到相关性 和聚类的索引
        maxXiangguanJL = max(minInC)
        maxXiangguanIndex = 0
        ii = 0
        while ii < len(julei):
            if minInC[ii] == maxXiangguanJL:
                maxXiangguanIndex = ii
                break
            ii = ii + 1
        #print('最相关的一个聚类：%d' % maxXiangguanIndex)
        #print('最大相似度：%f' % maxXiangguanJL)
        if maxXiangguanJL > menxian:
            # 属于已经存在的聚类  加入到该聚类
            julei[maxXiangguanIndex].append(int(i))
        else:
            # 新生成一个聚类
            temp = [i,]
            julei.append(temp)
    #print('聚类结果：')
    #print(julei)

    #模板选取原则：
    # 1、选出聚类数量前百分之三十的聚类
    # 2、分别计算前百分之三十聚类中 每条信号在len(sigModel/2)处的平均值  选择最大的那个聚类
    # 3、在选出的聚类中找出峰值最高的那条信号作为模板
    juleiLens = np.zeros(len(julei),dtype=np.int)
    i = 0
    while i < len(julei):
        juleiLens[i] = len(julei[i])
        i = i + 1
    num30 = int(len(julei)*0.3)
    maxNum = np.zeros(num30)
    maxNumIndex = np.zeros(num30,dtype=int)
    i = 0
    while i < num30:
        maxNum[i] = max(juleiLens)
        maxNumIndex[i] = np.argmax(juleiLens)
        juleiLens[maxNumIndex[i]] = 0
        i = i + 1
    # print('数量前百分之三十的聚类分别是：')
    # print(maxNumIndex)
    # print('数量分别是：')
    # print(maxNum)
    average = np.zeros(num30)
    i = 0
    while i < num30:
        ii = 0
        temp = 0
        while ii < maxNum[i]:
            temp = temp + beats[julei[maxNumIndex[i]][ii]][int(fs*0.4)]
            ii = ii + 1
        temp = temp/(maxNum[i])
        average[i] = temp
        i = i + 1
    maxIndex = np.argmax(average)
    maxjuleiIndex = maxNumIndex[maxIndex]
    # print('它们的平均值是：')
    # print(average)
    # print('平均值最大的聚类是：')
    # print(maxjuleiIndex)
    middleValeu = np.zeros(len(julei[maxjuleiIndex]))
    i = 0
    while i < len(middleValeu):
        middleValeu[i] = beats[julei[maxjuleiIndex][i]][int(fs*0.4)]
        i = i + 1
    sigModelIndex = np.argmax(middleValeu)
    sigModel = beats[julei[maxjuleiIndex][sigModelIndex]]
    # print('聚类中的中值分别是')
    # print(middleValeu)
    # print('最大的中值索引是：')
    # print(sigModelIndex)
    # print(sigModel)

    return julei,beats,sigModel

def computerCorrFunc(signalIn,sigModel):
    """
    计算相关系数函数  返回值长度等于signalIn长度
    :param signalIn: 输入信号
    :param sigModel: 信号模板  0.8秒
    :return: 相关系数函数
    """
    windowLen = len(sigModel)
    corrFunc = np.zeros(len(signalIn))
    i = int(windowLen/2)
    while i < len(signalIn)-windowLen/2:
        sig = np.zeros([2,windowLen])
        sig[0] = signalIn[int(i-windowLen/2):int(i+windowLen/2)]
        sig[1] = sigModel
        corrTemp = np.corrcoef(sig)
        corrFunc[i] = corrTemp[0,1]
        i = i + 1
    return corrFunc

def computerCorrFuncPecks(signalIn,sigModel):
    """
    计算相关系数函数  仅匹配极值点的相关系数  得到有关极值点的相关函数
    :param signalIn: 输入信号
    :param sigModel: 输出信号
    :return: 相关系数函数
    """
    windowLen = len(sigModel)
    peacks,_ = signal.find_peaks(signalIn)
    corrFunc = np.zeros(len(signalIn))
    i = 0
    for peack in peacks:
        if peack < windowLen/2 or peack > len(signalIn)-windowLen/2:
            continue
        else:
            sig = np.zeros([2,windowLen])
            sig[0] = signalIn[int(peack-windowLen/2):int(peack+windowLen/2)]
            sig[1] = sigModel
            corrTemp = np.corrcoef(sig)
            corrFunc[peack] = corrTemp[0, 1]
            i = i + 1
    return corrFunc

def findBeats(corrFunc,fs):
    """
    寻找心拍的开始点  算法流程：
    1、将corrFunc做三次方  加强信号区分度
    2、在前0.8秒 找到第一个值大于0.65的点  作为第一个心拍
    3、以上一个心拍点为起点  向后搜索0.5-1.2秒(125-300点)范围内的最大值点
    4、如果该点大于0  且心率变动不超过15% 则为有效心跳  否则选取次高点
    :param corrFunc:
    :param fs:
    :return:
    """
    beatsIndex = []
    corrFunc = corrFunc*corrFunc*corrFunc
    windowsLen = 1.2
    while True:
        winStart = int(fs * windowsLen) #在信号开始的1.2秒里寻找最值点
        signalTemp = corrFunc[0:winStart]
        maxIndex = np.argmax(signalTemp)
        if signalTemp[maxIndex] > 0.20:
            beatsIndex.append(int(maxIndex))
            break
        else:
            # 找不到有效心跳  重新匹配模板
            print('找不到有效心跳,增加搜索范围')
            windowsLen = windowsLen + 0.5
    i = 0
    while True:
        # print(len(beatsIndex))
        # print(beatsIndex[i])
        if len(corrFunc)-beatsIndex[i] > fs*0.5 :
            startIndex = int(beatsIndex[i] + fs*0.5)
            if len(corrFunc)-beatsIndex[i] > fs*1.2:
                endIndex = int(beatsIndex[i] + fs*1.2)
            else:
                endIndex = len(corrFunc)-1
            signalTemp = corrFunc[startIndex:endIndex]
            pecks, _ = signal.find_peaks(signalTemp)
            pecksValue = np.zeros(len(pecks))
            if len(pecks) < 1:
                break
            mm = 0
            while mm < len(pecks):
                pecksValue[mm] = signalTemp[pecks[mm]]
                mm = mm + 1
            if len(beatsIndex)>2:
                while True:
                    maxPecksValueIndex = np.argmax(pecksValue)
                    print(pecksValue)
                    if pecksValue[maxPecksValueIndex] < 0.01:
                        print('错误峰')
                        beatsIndex.append(min(int(startIndex + fs * 0.8),len(corrFunc)-1))
                        break
                    maxIndex = pecks[maxPecksValueIndex]
                    maxIndex = startIndex + maxIndex
                    diffRate = abs((maxIndex-beatsIndex[i])-(beatsIndex[i]-beatsIndex[i-1]))/(beatsIndex[i]-beatsIndex[i-1])
                    # print('maxIndex:%d' % maxIndex)
                    # print('beatsIndex[i]:%d' % beatsIndex[i])
                    # print('beatsIndex[i-1]:%d' % beatsIndex[i - 1])
                    # print('diffRate:%f' % diffRate)
                    if diffRate < 0.30:
                        beatsIndex.append(int(maxIndex))
                        break
                    else:
                        # 选取第二峰
                        pecksValue[maxPecksValueIndex] = 0
                        # print('选取第二峰')
                        # print('原始峰Index：%d' % maxIndex)
            else:
                # 选取最大的极大值点
                pecks2,_ = signal.find_peaks(signalTemp)
                jidaValue = np.zeros(len(pecks2))
                vv = 0
                while vv < len(pecks2):
                    jidaValue[vv] = signalTemp[pecks2[vv]]
                    vv = vv + 1
                if len(pecks2) < 1:
                    # print('错误峰')
                    beatsIndex.append(int(startIndex + fs * 0.8))
                jidaValueIndex = np.argmax(jidaValue)
                jidaIndex = pecks2[jidaValueIndex]
                beatsIndex.append(int(jidaIndex+startIndex))
                # print('pecks:')
                # print(pecks2)
                # print('jidaValue:')
                # print(jidaValue)
                # print('jidaValueIndex:')
                # print(jidaValueIndex)
                # print('jidaIndex')
                # print(jidaIndex)
                # print('index:%d' % int(jidaIndex+startIndex))
            i = i + 1
        else:
            break
    beatsIndex = np.array(beatsIndex)
    return beatsIndex

def findBreathBeats(signalIn , fs):
    """
    计算呼吸序列  使用简单的寻峰算法
    :param signalIn: 输入信号  经过带通滤波[0.1-0.8Hz]之后的信号
    :param fs: 采样率
    :return: 呼吸峰值点的索引序列
    """
    # 算法步骤：
    # 1、//输入信号归一化
    # 2、计算所有的极大值点
    # 3、在0-5秒的范围内找出第一个幅度超过0.5的峰  如果找不到则使用最高峰代替
    # 4、向后在1.25秒-6秒内按照峰的大小循环  幅值大于前一个值的60%即为有效点
    bcgmax, bcgmin = signalIn.max(), signalIn.min()
    signalIn = (signalIn - bcgmin)/(bcgmax-bcgmin)
    breathPause = 0
    limit = 10
    signalLength = len(signalIn)
    pecksRel = []
    pecks,_ = signal.find_peaks(signalIn)
    if signalLength < 5*fs:
        return pecks
    signalTemp = signalIn[0:5*fs]
    maxIndex = np.argmax(signalTemp)
    pecksRel.append(int(maxIndex))
    i = 0
    windowsLen = 5
    while True:
        if int(pecksRel[i]+fs*1.2) < signalLength:
            startIndex = pecksRel[i] + fs*1.2
            if int(startIndex+fs*windowsLen) < signalLength:
                endIndex = startIndex + fs*windowsLen
            else:
                endIndex = signalLength - 1
            # print('<------------------------------->')
            # print('本次计算区间：[%d:%d]' % (int(startIndex), int(endIndex)))
            # print('windowsLen:%d' % int(windowsLen))
            signalTemp = np.zeros(int(endIndex-startIndex))
            hh = 0
            while hh < len(signalTemp):
                signalTemp[hh] = signalIn[int(startIndex+hh)]
                hh = hh + 1
            #signalTemp = signalIn[int(startIndex):int(endIndex)]
            zeroIndex = 0
            hh = 0
            while hh < len(signalTemp):
                if signalTemp[hh] < 0.4:
                    zeroIndex = hh
                    break
                else:
                    hh = hh + 1
            hh = 0
            print('zeroIndex:%d' % zeroIndex)
            while hh < zeroIndex:
                signalTemp[hh] = 0.5
                hh = hh + 1
            pecks2,_ = signal.find_peaks(signalTemp)
            pecksValue = np.zeros(len(pecks2))
            ii = 0
            while ii < len(pecks2):
                pecksValue[ii] = signalTemp[pecks2[ii]]
                ii = ii + 1
            cc = 0
            print('本区间峰值索引及数值：')
            print(pecks2+startIndex)
            print(pecksValue)
            if windowsLen == 10:
                breathPause = breathPause + 1
                print('呼吸暂停：%d' % int(breathPause))
                limit = 50
            if windowsLen == 50:
                break
            if endIndex == (signalLength-1):
                break
            if len(pecks2) == 0:
                windowsLen = windowsLen + 1
                continue
            if windowsLen < limit:
                maxIndex = np.argmax(pecksValue)
                print('本区间最大值的索引为：%d' % int(pecks2[maxIndex]+startIndex))
                val = 0
                nn = 0
                while nn < 5:
                    if (i - nn) < 0:
                        break
                    val = val + signalIn[int(pecksRel[i - nn])]
                    nn = nn + 1
                aveVal = val/nn
                diff = abs(aveVal - pecksValue[maxIndex])/signalIn[int(pecksRel[i])]
                print('diff:%s' % str(diff))
                if diff < 0.8:
                    pecksRel.append(int(pecks2[maxIndex]+startIndex))
                    print('找到呼吸点：%d' % int(pecks2[maxIndex]+startIndex))
                    windowsLen = 5
                    limit = 10
                    i = i + 1
                else:
                    windowsLen = windowsLen + 1
        else:
            break
    return np.array(pecksRel)

def CorrectHeartBeatError(filtedHeart,beats,fs):
    """
    矫正心拍的误差  根据beats序列反向寻找原信号filtedheart中的峰值点  并返回新的pecks
    搜索范围：正负20ms
    :param filtedHeart:滤波后的心拍信号
    :param beats:原始序列
    :param fs:采样率
    :return:新的新拍序列
    """
    i = 0
    pecks = []
    while i < len(beats):
        if beats[i]-fs*0.04 > 0:
            startIndex = int(beats[i]-fs*0.04)
        else:
            startIndex = 0
        if beats[i]+fs*0.04 < len(filtedHeart)-1:
            endIndex = int(beats[i]+fs*0.04)
            signalTemp = filtedHeart[startIndex:endIndex]
            maxIndex = np.argmax(signalTemp)
            pecks.append(int(startIndex+maxIndex))
            i = i + 1
        else:
            break
    return np.array(pecks)


def computerHeartBeat2Hz(HeartBeats,fs,signalLength):
    """
    对输入的心拍信号进行三次样条插值  以2Hz的频率重采样
    :param HeartBeats: 心拍信号  每个J峰的索引
    :param fs: 输入信号的采样率
    :return: 2Hz重采样的信号  长度为signalLength*2
    """
    HeartBeats0 = np.zeros(len(HeartBeats)+1)
    i = 0
    while i < len(HeartBeats):
        HeartBeats0[i+1] = HeartBeats[i]
        i = i + 1
    HeartBeatsy = np.diff(HeartBeats0)*(1000/fs)
    HeartBeatsy[0] = HeartBeatsy[1]
    f1 = interpolate.splrep(HeartBeats,HeartBeatsy)
    x = np.linspace(0,signalLength*fs,signalLength*fs)
    y1 = interpolate.splev(x,f1)

    i = 0
    y2 = np.zeros(signalLength*2)
    while i < signalLength*2:
        y2[i] = y1[int(i*(fs/2))]
        i = i + 1

    return y2

def findBreathBeatsNew(signalIn ,fs):
    signalLength = len(signalIn)
    pingjunLine = np.zeros(signalLength)
    if signalLength < fs*30:
        return np.zeros(1)
    signalTemp = np.zeros(150)
    i = 0
    while i < 150:
        signalTemp[i] = signalIn[i*50]
        i = i + 1
    signalpinpu = np.fft.fft(signalTemp)
    pinlv = np.argmax(signalpinpu[0:70])*0.0333
    if pinlv == 0:
        pinlv = 0.3
    T = int((1/pinlv)*250)
    i = 0
    while i < T:
        pingjunLine[i] = np.average(signalIn[0:T*2])
        i = i + 1
    i = 0
    while i < signalLength-2*T:
        pingjunLine[i+T] = np.average(signalIn[i+T:i+2*T])
        i = i + 1
    i = 0
    while i < T:
        pingjunLine[signalLength-T+i] = np.average(signalIn[signalLength-2*T:signalLength])
        i = i + 1
    signalfu = np.zeros(signalLength)
    i = 0
    while i < signalLength:
        signalfu[i] = -(signalIn[i]-pingjunLine[i])
        if signalfu[i] < 0:
            signalfu[i] = 0
        i = i + 1
    pecks = []
    i = 0
    startIndex = 0
    endIndex = 1
    while i < signalLength:
        if signalfu[i] == 0:
            endIndex = i
            signalTemp1 = signalfu[startIndex:endIndex]
            if len(signalTemp1) == 0:
                i=i+1
                continue
            pecks.append(np.argmax(signalTemp1)+startIndex)
            # print(pecks)
            ii = 0
            while ii < signalLength-endIndex:
                if signalfu[endIndex+ii]>0:
                    startIndex = endIndex + ii
                    i = startIndex
                    break
                else:
                    ii = ii +1
            if ii == signalLength-endIndex:
                break
        else:
            i = i + 1
    pecks = np.array(pecks)
    # plt.plot(signalfu)
    # plt.plot(pecks,signalfu[pecks],'*')
    # plt.show()



    return pingjunLine,pecks
